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Abstract 

Plasticity is governed by the evolution of, in general anisotropic, systems of dislocations. We seek to faithfully 
represent this evolution in terms of density-like variables which average over the discrete dislocation microstructure. 
Starting from T. Hochrainer’s continuum theory of dislocations (ODD) [Hochrainer 2015], we introduce a methodology 
based on the ’Maximum Information Entropy Principle’ (MIEP) for deriving closed-form evolution equations for 
dislocation density measures of different order. These equations provide an optimum representation of the kinematic 
properties of systems of curved and connected dislocation lines with the information contained in a given set of density 
measures. The performance of the derived equations is benchmarked against other models proposed in the literature, 
using discrete dislocation dynamics simulations as a reference. As a benchmark problem we study dislocations moving 
in a highly heterogeneous, persistent slip-band-like geometry. We demonstrate that excellent agreement with discrete 
simulations can be obtained in terms of a very small number of averaged dislocation fields containing information 
about the edge and screw components of the total and excess (geometrically necessary) dislocation densities. From 
these the full dislocation orientation distribution which emerges as dislocations move through a channel-wall structure 
can be faithfully reconstructed. 

Keywords: continuum theory of dislocations, dislocation dynamics, persistent slip bands, 
alignment tensors 


1. Introduction 

The development of physically based and predictive models for plasticity on the micro meter 
scale has been a challenging task ever since the connection between dislocations and the macroscopic 
material response has been made in the 1930s. On the one hand, continuum approaches based on 
naive averaging for obtaining continuous densities without taking into account the curved and line¬ 
like nature of systems of dislocations must fail in many cases because important microstructural 
information is being lost; discrete approaches, on the other hand, contain full information about 
the microstructure. However, the number of dislocations or accumulated plastic strain can be high 
which can easily become a limiting factor even for today’s discrete dislocation dynamics (DDD) 
models. Continuum models of dislocation systems which average over the discrete dislocation 
microstructure do not suffer from this restriction and hence might have advantages over discrete 
models, provided the motion of dislocations can be correctly captured in such a continuum (or 
dislocation-density based) framework. 

Historically, dislocation density-based plasticity models have evolved along several independent 
lines. The first line originates from the continuum theory of dislocations and internal stresses as 
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developed by Kroner (19581 and Nye (19531. This theory is formulated in a geometrically rigor¬ 
ous manner and provides generic relationships between the dislocation microstructure, the plastic 
distortion and the associated internal stress fields. The fundamental object of the theory is the 
dislocation density tensor a which is defined as the curl of the plastic distortion, a = —curl^SP^ 
This theory was extended by Mura (1963) who formulated a kinematic equation of evolution for 
the dislocation density tensor, dtcx = —curl[r; x q] where v is the dislocation velocity vector. In this 
form the theory provides a full description of dislocation microstructure evolution and of plastic 
deformation for situations where all dislocations are geometrically necessary dislocations (GND), 
i.e., where they can be envisaged as contour lines of the plastic shear strain on the respective slip 


systems (e.g. Sedlacek et al., 2003; Xiang 2009 Zhu and Xiang, 2015). In the general case where 


dislocations of multiple orientations and slip systems are present, such a field theory, in order to 
fully capture the evolution of the dislocation microstructure, requires a spatial resolution that is 
well below the spacing of the individual dislocation lines in order to make the dislocation velocity 
field V uniquely defined. (See e.g. Xia and El-Azab (2015) and Zhang et al. (2015) for such imple¬ 
mentations). A similar spatial resolution is also required for phase field approaches to dislocation 


microstructure evolution who directly simulate the evolution of the shear strain fields (Wang et al. 


2001 Rodney et al., 2003). Such simulations, which evidently incur a high computational cost. 


may be envisaged as field-theoretical approaches to discrete dislocation dynamics simulation and 
will not be discussed in the following. Instead our focus of interest is on average, statistical de¬ 
scriptions of the dislocation microstructure which work in situations where dislocations of multiple 
orientations are present within the same volume element, such that a mapping on corresponding 
density-like variables provides an efficient compression of information. On such a coarse grained 
level, approaches which are directly based on spatial averaging of the classical dislocation density 
tensor a can work only in exceptional circumstances: Once the fine structure of the plastic strain 
field is smoothed over, the dislocation lines associated with the averaged out features are no longer 
represented by the dislocation density tensor. Unfortunately, this background of ’’statistically 
stored” dislocations still contributes both to plastic flow and to work hardening. We are thus faced 
with the fundamental problem how they can be incorporated into a density-based theory. 

The second line of continuum models, which originates from the work of [Johnston and Gilman 


(1959) and Kocks (1976), solved this problem by taking the radical approach of disposing with 
geometry altogether and considering only the ’’statistically stored” contribution to the dislocation 
density, which is characterized by a scalar density measure p of dimension 1/length^. For this 
density measure, phenomenological evolution equations are formulated which relate the dislocation 
density to the strain, p = ^( 7 ), and these equations were combined with other phenomenolog¬ 
ical relations which relate the scalar density of dislocations to the flow stress (e.g., the Taylor 
relationship r = aGhp^/"^ where G is the shear modulus, b the Burgers vector modulus, and a a 
dimensionless constant) and to the plastic flow rate (e.g. the Orowan relation dj/dt = pbv where 
V is the scalar magnitude of the dislocation velocity). This modelling approach and its derivatives 
were successfully used to formulate phenomenological models of work hardening and plastic flow 
(e.g. Estrin and Mecking 1984 Gaceres and Blake, 2007 Bouaziz et al. 2013). Generally speaking. 


the approach works well as long as deformation is, on the scale of the description, homogeneous 
such that geometrically necessary dislocations - or, equivalently, strain gradients - need not to be 
taken into accounlQ 

In recent years, several attempts have been made to unify both strands of continuum modeling 


^The limitations of the line of Gilman and Kocks can easily be seen: consider, for instance, the deformation of 
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and to arrive at models which can capture the combined evolution of ’’statistically stored” and 
’’geometrically necessary” dislocation densities in a framework which can faithfully represent the 
underlying motion of the discrete dislocation lines. Pioneering work was done on two-dimensional 


(2D) systems of straight parallel dislocations by Groma (1997|; Zaiser et al. (20011; Groma et al. 


( 2003[ ) who systematically formulated evolution equations for 2D dislocation densities as statistical 
averages over the dynamics of the corresponding discrete dislocation systems. Generalizing this 
approach to three-dimensional systems of curved dislocations has, however, proven to be quite 
challenging: Straight parallel dislocations can be envisaged as points in the intersecting plane, and 
both the mathematical definition of densities of such objects and a consistent formulation of their 
kinematics are almost trivial tasks. Three-dimensionally moving dislocations, on the other hand, 
are curved lines which move perpendicular to their line direction while remaining topologically 
connected, and the definition of densities of such objects and corresponding formulation of their 
kinematics is far from straightforward. As a consequence, models have been (and up to date still 
are) developed that intentionally neglect some aspects: e.g. screw-edge representations (Arsenlis 


et al., 2004 Reuber et al., 2014 Leung et al. 20151 are a coarse way of approximating continuously 


curved dislocation loops, while other approaches rely on a geometrically consistent description of 
the evolution of the GND content as described by the Kroner-Nye density tensor but complement 
this with phenomenological ad-hoc assumptions regarding the evolution of the statistically stored 
dislocation density and/or the contribution of such dislocations to the plastic strain rate (e.g. 
Acharya and Roy, 2006[ [Fressengeas et aT 20111. 


We finally note another line of continuum dislocation theories which is based on the low-energy 
dislocation structure (LEDS) hypothesis proposed by Hansen and Kuhlmann-Wilsdorf (19861. The 
LEDS-hypothesis to a certain extent dismisses the role of the history (initial dislocation microstruc¬ 
ture and deformation path) and indicates that among all admissible dislocation configurations the 
dislocation microstructure emerging under a given set of boundary constraints is the one that min¬ 
imizes the energy of the crystal. Along this line Le and Gunther (2014) developed a continuum 
dislocation theory (GDT) of straight dislocation lines. Recently Le (2016) generalized this theory 
to 3D curved dislocation lines (3D-GDT) using scalar density like variables that contain informa¬ 
tion about the dislocation line orientation and curvature - an approach that has strong analogies 
to our higher-dimensional continuum theory discussed below. The 3D-GDT has shown promise 
in reproducing analytical solutions of test cases such as dislocation pile ups in simple shear de¬ 
formation. However the theory is formulated under the assumption that “the dislocation network 
is regular in the sense that nearby dislocations have nearly the same direction and orientation” 


(Le, 2016) which implies that the dislocation microstructure should locally consist of GNDs only. 
This assumption restricts the applicability of the theory - in particular, it is difficult to see how it 
could be applied to situations where each averaging volume contains exclusively or predominantly 
dislocations of zero net Burgers vector, or more generally speaking to the interplay of SSD and 
GND evolution which is the main focus of the present work. 

Eirst important steps towards representing generic 3D systems of curved and connected dislo¬ 


cation lines through density measures date back to the work of Kosevich in the 1970s (Kosevich 


a polycrystal before grain boundaries yield plastically. Any approach which considers the dislocation density an 
increasing function of the strain will predict that dislocation densities are largest in the grain interior, where the 
strains have their maximum. However, it is obvious that in reality dislocation densities will be largest near the grain 
boundaries where dislocation motion is constrained and dislocations form pile ups which accommodate the associated 
strain gradients. 
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19791 and were later taken up by El-Azab (20001. These authors introduced the idea of envisaging 


dislocations in a phase or state space where densities carry additional information about their line 
orientation. This approach was systematized and formulated in a mathematically rigorous manner 
by Hochrainer and co-workers who defined dislocation densities and derived their evolution in a 
higher-dimensional configuration space containing line orientation variables as extra dimensions 
(Hochrainer, 2006; Hochrainer et alT| 2007[ Sandfeld et al. 2010 2015a). This higher-dimensional 
continuum dislocation dynamics (hdCDD) theory is based upon statistical averaging of the kine¬ 
matics of dislocations, i.e., it provides a density based representation of the manner how curved 
and connected lines move in a given velocity field, with results that have been demonstrated to be 


in quantitative agreements with those derived from the corresponding discrete dynamics (Sandfeld 


et al. 2015a I. The theory can easily handle anisotropic orientation distributions and direction 


dependent dislocation velocities (Sandfeld et al., 2015b|. The drawback is the high computational 


cost of implementing a dynamics on a higher-dimensional space where, in each spatial volume ele¬ 
ment, one needs to consider the evolution of the continuous orientation distribution of dislocations. 

As a consequence, attempts have been made to formulate simplified variants of the theory (all 
denoted as CDD) (Hochrainer et al. 2009 Sandfeld et al. 2011 Hochrainer et al. 2014 Hochrainer[ 
20151 which consider density-like variables that represent low-order moments of the dislocation ori¬ 


entation distribution only. Such CDD models were applied to a number of benchmark problems 


Hochrainer 


f size effects (Sandfeld et al. 

2011 Sandfeld and 

Sandfeld and Zaiser 

2015). 

Hochrainer 

(2015) 


put these attempts on a systematic foundation by observing that the information contained in 
the orientation-dependent dislocation density function is exactly the same as that contained in an 
alignment tensor expansion of this function with respect to the orientation variables. The com¬ 
ponents of the dislocation density alignment tensors can be envisaged as density-like fields which 
contain more and more detailed information about the orientation distribution of dislocations. For 


the dislocation density alignment tensors, Hochrainer (2015) derived an infinite hierarchy of evo¬ 


lution equations where the evolution of alignment tensors of low order depends on higher-order 
alignment tensors. In the present paper we present a systematic method for closing this hierarchy 
at any given order, and we assess the performance of the resulting lowest-order CDD models. 

The goal of this paper is threefold: (i) to introduee a systematie approaeh for elosing the evolu¬ 
tion equations of CDD using the Maximum Information Entropy Prineiple (MIEP); (ii) to ealeulate 
the elosure approximations for the two lowest order variants of CDD and diseuss the eonsequenees 
of these approximations; (Hi) to benehmark the performanee of the resulting CDD models in terms 
of their ability to eorreetly represent the way how disloeations move in a strongly anisotropie dislo- 
eation mierostrueture. 


In Section we start by introducing the notions of kinematic and dynamic consistency as 
measures for the performance of averaged density-based theories of dislocation systems. This is 


followed by a brief overview of the hdCDD theory in Section 3.2 In Section 3.3 generic steps 


(based on alignment tensors) for deriving simplified CDD models with ’customize-able’ accuracy 
from hdCDD are outlined. We then turn in Section 3.4 to the work horse of this paper, the 
Maximum Information Entropy Prineiple (MIEP), which will be used both as a tool for deriving 
CDD evolution equations, and as a means of analyzing what information is lost in simplified CDD 
versions as compared to hdCDD. In Sections |3.5| - |3T| we then derive the two simplest possible CDD 
variants. Their performance in terms of the so-called ’kinematic consistency’ (i.e. how well the 
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evolution of curved and connected dislocations in a given velocity field is represented) is numerically 
benchmarked and analyzed in Section]^ where we compare them to results from DDD and hdCDD 
simulations, as well as to some of alternative models published in the literature. As a benchmark 
example we use a problem that has been extensively studied also in the experimental literature, 
namely the motion of dislocations in the periodic wall structure of persistent slip bands in fatigued 
metals. 

2. Kinematic and dynamic consistency: How to measure the quality of averaged de¬ 
scriptions of dislocation systems 

The first essential step in constructing an averaged continuum theory of dislocations - even 
though this is rarely made explicit - consists in the definition of statistical rules which map a 
discrete dislocation system (or an ensemble of such systems) onto a set of continuous density 
measures. Furthermore, we need to formulate equations of evolution for the density measures which 
map an initial set of density measures onto a set of measures at a later time. The performance of 
a continuum theory may then be assessed by investigating to which extent this mapping meets the 
requirements of kinematic and dynamic consistency with the underlying discrete dynamics: 

The kinematic consistency of a continuum theory of dislocation motion can be envisaged as 
its ability to correctly account for the fact that dislocations are curved and connected lines which 
evolve by moving perpendicular to their line direction. It is thus a measure of the ability of a 
continuum theory to still account, after averaging, for the peculiar geometrical properties of the 
underlying discrete objects. To assess kinematic consistency, we propose to use a benchmarking 
method: For benchmark problems, we construct an initial discrete configuration (or an ensemble 
of such configurations) and evaluate the corresponding density measures. We then impose a dislo¬ 
cation velocity field v{r) where the magnitude of the dislocation velocity does not depend on the 
dislocation configuration - i.e., we neglect internal stresses and dislocation interactions (obtaining 
those would be the problem of dislocation dynamics). We investigate in parallel the evolution of 
the discrete system by moving all lines with velocity v{r) perpendicular to their direction, and the 
evolution of the density measures under the respective evolution equations. Finally, after some 
time t, we evaluate the density pattern resulting from the evolved discrete system and compare 
with the evolved continuous measure. Kinematic consistency of a continuum theory can then be as¬ 
sessed in terms of the degree of agreement between the evolved continuous system, and the density 
representation of the evolved discrete system (c.f. Sandfeld and Po (20151). 

The dynamic consistency of a continuum theory concerns its ability to predict the actual evolu¬ 
tion of dislocation systems under external load. Evaluation proceeds in much the same manner as 
for dynamic consistency, however, the dislocation velocity field is not imposed externally. Instead, 
a set of external boundary conditions is specified and the dislocations evolve under the influence 
of the associated stresses, as well as of their mutual interactions. The local dislocation velocities 
are computed as functionals of stress and of the variables describing the dislocation configura¬ 
tion. Consistency is again evaluated by comparing the evolved continuous system with the density 
representation of the evolved discrete system(s). 

Kinematic consistency may be considered a special case of the more general concept of dynamic 
consistency. The two concepts coincide in certain limiting cases where the dislocation velocity is 
mainly controlled by ’external’ stresses while dislocation interactions are unimportant - an example 
being semiconductor single crystals with ultra-low dislocation densities (and thus small dislocation 
interaction stresses) and high Peierls barriers. In such systems we may write the dislocation 


Sandfeld and Po (2015 
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velocities as functions of the stresses arising from external boundary constraints or differential 
thermal expansion, without consideration of variables characterizing the dislocation microstructure. 
In general, however, dynamic consistency is a much more stringent and demanding requirement 
than kinematic consistency. While kinematic consistency merely tells us that a continuum model 
accounts for the geometrical nature of dislocations, dynamic consistency additionally requires that 
the model correctly represents the physics of dislocation interactions, i.e. the manner how internal 
stresses and the resulting long-range and short-range dislocation interactions govern dislocation 
motion. So the reader might legitimately ask why we are so obsessed with kinematics. 

There are two answers to this question. First, we insist that kinematic consistency is a neces¬ 
sary prerequisite for achieving dynamic consistency: Without investigating the issue of kinematic 
consistency we do not even know whether an averaged theory correctly captures where dislocation 
lines are moving (kinematic consistency), which makes it quite futile to ask the question how fast 
they are moving (dynamic consistency). We note that, in order to ensure kinematic consistency, 
it is not sufficient to ensure that the dislocation density tensor a, as calculated from the density 
measures, is divergence free (divo: = 0). This requirement, which ensures that the geometrically 
necessary dislocations can be considered a set of closed contour lines of the strain field, is a nec¬ 
essary condition for kinematic consistency but not sufficient to ensure correct representation of 
the coupled kinematics of systems containing both statistically stored and geometrically necessary 
dislocations. Second, and possibly more important, solving the problem of dislocation kinematics 
is a step which can help to systematically derive theories of dislocation dynamics. In a recent paper 
of Hochrainer (20161 it was demonstrated that, if the energy of a dislocation system is known as 
a functional of a set of density measures and the kinetic equations which govern the evolution of 
the same measures are also known, then one can use the powerful framework of linear irreversible 
thermodynamics to systematically construct thermodynamically consistent dynamic theories. In 
another recent paper by Zaiser (2015) it was further demonstrated that the elastic energy functional 
of a discrete dislocation system can be averaged to express the elastic energy as a functional of the 
corresponding dislocation density alignment tensors - that is, of the CDD field variables formu¬ 
lated by Hochrainer (20151 and also used in the present work. Knowledge of the averaged energy 
functional together with a correct average representation of the kinematics means that a consistent 
dynamic theory can be derived without any need for making arbitrary constitutive assumptions. 
The present paper aims at finalizing the kinematic side of this project. 


3. Theory 

3.1. Mathematical notations and conventions 

In the following x denotes the cross product between two vectors, a dot denotes a single 
contraction and ® is the tensor product. We consider dislocations moving on a single slip system 
with Burgers vector b and slip plane normal vector n. Calculations are performed using a Cartesian 
coordinate system with basis vectors 61 , 62,63 where without loss of generality we set 61 = b/\b\, 
62 = (i/\a\ with a = n X b and 63 = n. We consider dislocations moving by glide only, hence 
the dislocation lines are contained in parallel glide planes. The spatial direction of the dislocation 
line can thus be characterized by a single angle (p, which we take to be the angle between the line 
tangent vector I and the Burgers vector such that l((p) = cos((^) 6 i -|-sin((/j) 62 . If l{(p) points into 
the direction of motion as we move counter-clockwise around a dislocation loop, the loop is termed 
positive, otherwise it is termed negative. By convention, positive loops expand and negative loops 
shrink under positive resolved shear stress. 
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The operator O'*- defined by t-*- = t x n rotates a vector t that is contained in the slip plane by 
90° around the slip plane normal. This operator can alternatively be expressed as t-^ = t.e where 
e is the second order Levi-Civita tensor in the slip plane coordinate system, with components = 
where are the components of n (we use the Einstein summation convention throughout). 
A partial derivative w.r.t. Xi is denoted by di = The nabla operator of spatial derivatives is 
defined as the vector operator V = Sidi, the curl operator is curl = Vx. 

^th completely symmetric tensor obtained from the vector t according to the 

recursion relation = t, (g) t. is the order symmetric alignment tensor of 

an orientation distribution function and defined as = j> ^((p)l((f)^"^d(p . Tr(T) gives 

the trace of a symmetric tensor T through summation over any two indices of T. The operator 
[T]sym completely symmetrizes a tensor by averaging its components over all possible permutations 
of indices. 


3.2. The higher dimensional eontinuum theory of disloeations 

Here, we only present a brief summary of the higher dimensional continuum theory of disloca¬ 


tions; for details the reader is referred to the published literature (Hochrainer et al. 2007 Sandfeld 


et al. 2010 Hochrainer et al.[ 20141. Dislocation lines are envisaged as lifted curves in a configu¬ 
ration space where each point {r,ip) contains the spatial point r and the orientation angle (p. To 
describe the kinematics of single lines within this space we introduce a generalized line direction 
L and a generalized velocity V, which denote the tangent to the lifted line and the velocity of the 
lifted line, respectively: 


L{r,ip) = {cosip,smip,k{r,ip)) = {l{ip),k{r,ip)) (1) 

V (r, if) = (usin (^, —ucos ip, i?(r, cp)). (2) 

Here I is the spatial line direction, v is the modulus of the spatial velocity v (the first two terms 
in (|^) which is normal to I, and k is the line curvature. The additional component of the lifted 
velocity vector, 'd{r,(p) = —Viu(r,(/j), is a rotational velocity which causes a line to move in 
ip direction, i.e. to change its orientation. denotes the gradient along the generalized line 
direction. With the above notations it is possible to define the so-called ’dislocation density tensor 
of second order’ along with its evolution equations. This tensor can be expressed in terms of a 
density function p{r, ip) which gives the density of dislocations for each line orientation separately 
as 

OL^{r,ip) = p{r,ip)L{r,ip)®h. (3) 

The generalized line direction L contains, besides the spatial line direction I, as additional compo¬ 
nent the local curvature k. Accordingly, the corresponding components of the higher order dislo¬ 
cation density tensor contain a new variable, the product of the orientation-dependent dislocation 
density and the local curvature, which we denote as curvature density g(r, ip) = p{r, ip) k{r, ip). In 
analogy to the condition that a is divergence free, one can show that also o:^ is solenoidal, i.e., 

Div = 0 44- cos ipdx{p) -h sin ipdy{p) -h dy>{q) = 0, (4) 

where Div(») := dx{») + dy{») + 5<^(»). (Q reflects the physical fact that dislocation lines do not 
start or end inside a crystal. The evolution equation for is 

dta^{r,ip) = —c\iT\{V{r,ip) x Q^{r,ip)). (5) 
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This evolution equation for Q^{r, (p) can be re-written in terms of two evolution equations for the 
scalar dislocation density p{r,(p) and the curvature density q{r,ip): 

dtp = -Div(pF) qv (6) 

dtq = — 0 ^( 9 ^ — pLi}) (7) 

In Q, the first term governs the transport of scalar density in the configuration space, while the 
last term is a source term and accounts for the change of density due to the expansion or shrinkage 
of loops. Similarly, the terms in Q govern the transport of the scalar curvature density in the 
configuration space. 


3.3. Reducing the higher dimensional continuum theory 

If all dislocations contained in a volume element share the same orientation, and thus move in 
the same direction, the extension to a higher dimensional configuration space which is at the core 
of hdCDD is redundant. However, the advantages of Hochrainer’s formulation become manifest 
as soon as we consider averaged dislocation fields where dislocations of different orientation are 
contained within the same volume element. Since the differently oriented segments inhabit differ¬ 
ent locations in configuration space, their motion in different directions can be resolved without 
problems. However, this capability comes at a very substantial price: First, the set of evolution 
equations, (§ and Q are from a numerical point of view difficult so solve due to the coupling of 
spatial and orientation degrees of freedom. Second, the additional orientation dimension introduces 
a large number of additional degrees of freedom, because this direction needs to be discretised such 
that gradients along the orientation direction can be properly resolved. This led to efforts towards 
simplifying the higher-dimensional evolution equations while preserving their performance in cap¬ 


turing dislocation kinematics (Hochrainer et al. 2009 Sandfeld et al. 2011 Hochrainer et ah 


20141. Such simplifications can be constructed in a systematic manner by performing a Fourier 


expansion of the functions p{r,(p) and q{r,(f) and using the resulting Fourier coefficients as vari¬ 


ables of a simplified continuum theory (CDD). Recently Hochrainer (20151 revisited this problem 


and introduced a mathematically more concise and elegant formulation which uses an expansion 
of the dislocation density and curvature density function in terms of symmetric alignment tensors. 
In the following, we just state the final results. Assuming pure glide motion for dislocation lines, 
the reducible symmetric alignment tensors of the dislocation and curvature densities are defined 
as: 


pW(r) j^p{r,p)l{ip)^'^ dip. 

= j>q{r,p) [e ■ l{p) ® dp. 


( 8 ) 

( 9 ) 


The symmetric alignment tensors are density-like quantities which characterize the dislocation 
microstructure. The zeroth order alignment tensors of p{r,p) and q{r,p) are: 


/9*(r) := p^^\r) := j>p{r, p) dp and q^{r) := q^^\r) := j>q{r, p) dp. 


( 10 ) 


Their physical meaning is as follows: p* gives total dislocation line length per unit volume, irre¬ 
spective of orientation. The zeroth-order alignment tensor q^ is the line length per unit volume, 
divided by the line curvature radius. If averaged over a volume of sufficient size, this quantity can 
alternatively be envisaged as a volume density of dislocation loops. 
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The first order dislocation density alignment tensor is a vector which measures the excess 
(geometrically necessary) dislocation density : 

p^^\r) = j^p{r,ip)l{ip)d(p ( 11 ) 

The components of this dislocation density vector p^^\r) = [/3i^^(^’),/ 02 ^^(?’)] are the screw and 
edge components of the GND density. p^^Velates to the curl of the plastic distortion (arising from 
slip on the single considered slip system), and hence to the classical dislocation density tensor, by: 

(g) 6 = Q = —curl^QP* (12) 

The second-order dislocation density alignment tensor is given by 

p^‘^\r) = jp{r,ip)l{ip)®l{ip)dp. (13) 

( 2 ) ( 2 ) 

The components p)/ and P 22 of this tensor can be understood as total densities of screw and edge 
components of the dislocation lines. Their sum, Tr(p(^^), is equal to p*. This is an example of the 
general relation 

Trp(") = and Trg(”) = (14) 

n 

Higher-order dislocation density alignment tensors are less straightforward to interpret in phys¬ 
ical terms. However, we think that the above defined quantities provide a sufficiently accurate 
statistical characterization of dislocation microstructures that accounts for all physically relevant 
distinctions: The distinction between screw and edge orientations (relevant because of potentially 
different physical properties such as mobility and cross slip) and the distinction between disloca¬ 
tions associated with strain gradients (the dislocation density vector which relates to the curl of the 
plastic distortion, or equivalently to the classical dislocation density tensor) and dislocations which 
are not. In particular, a theory using only the above quantities is perfectly capable of handling 
situations where locally the dislocation density tensor (the averaged curl of the plastic distortion) 
is zero but the total, screw and edge dislocation densities are not. 

The evolution equations for p^"^) and form an infinite hierarchy where the evolution of 
lower-order tensors depends on higher order tensors. In the case of an isotropic dislocation velocity 


they are given by 

dtp^^^ = V • {ve ■ p*-^^) -|- (15) 

dtp(^^ = [v • (ve 0 p("-i)) + (n - - (n - l)e ■ ■ Vv] (16) 

5^q(0) = V • (ug(^) - p(2) . Vv) (17) 

= [V • (ug(”+^) - £ • (18) 

L J sym 

where are auxiliary symmetric curvature tensors defined as 

g("-)(r) j)q(r, ip) e ■ l{(p) (g) e ■ l{ip) (g) (ig^ 
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In order to close the hierarchy of evolution equations for the dislocation density alignment tensors 
at order n we need to express and in terms of lower order tensors. We note that a 

matching closure relation for the curvature density alignment tensors is not required, because 
for n > 1 they can be derived from those for p using the solenoidality condition 0. This yields 


^W(r) = - [divp(’^+^)(r) 
n L 


Jsym 


( 20 ) 


The only variable which cannot be expressed in this manner is the zeroth-order alignment tensor 
q(0) = which must be considered as an independent dynamic variable of the theory. 


3.4- Recovering the dislocation orientation distribution function (DODF) based on the Maximum 
Information Entropy Principle (MIEP) 

In Section |3.3| we concluded that the evolution equations of the dislocation density function 
p(r, ip) can be expressed in terms of its alignment tensors. The information about dislocation 
orientations is contained in the “dislocation orientation distribution function” (DODF) which is 
defined as Pr(<p) = p)Ip^ and gives the probability density for dislocations to have, in a given 
spatial point, a certain orientation ip. In the following we shall drop, for simplicity of notation, the 
explicit dependency on r with the understanding that p must be considered in each spatial point 
separately. 

The infinite set of alignment tensors contains complete information about the orientation dis¬ 
tribution function, and conversely, closing the hierarchy at any given order n results in a loss of 
information: Instead of a single orientation distribution, we now have an infinite number of pos¬ 
sible distribution functions which are consistent with the known low-order alignment tensors. To 
estimate the full distribution function based on the limited information contained in the known 
alignment tensors, one may use the Maximum Information Entropy Principle (MIEP) which states 
that among all possible distributions the most probable is the one that maximizes the associated 
Shannon information entropy defined as: 

H=^ - j> p{ip) ln{p{ip))dip, (21) 


The ’’best guess” for p{ip) obtained by maximizing H needs to be consistent with the information 
contained in the n known dislocations alignment tensors, hence it must fulfill the constraints 


j>p{ip)dip = l, ... , 


d^ = pW/ph 


( 22 ) 


These constraints are however not independent since, from the alignment tensors of order n and 
n — 1, all lower-order alignment tensors can be constructed according to Eq. (14). Thus, we 
are only dealing with 2n -|- 1 independent constraints. Eor these we could consider the 2n -|- 1 
independent components of the completely symmetric tensors and However, for our 

later considerations it is more convenient to use an alternative choice and define 2n-|-1 independent 
constraints in terms of the first two components p^ii i and p^ 2 i i the alignment tensors at each 
order k < n: 


a) n -|- 1 constraints, for all 0 < A: < n : 

b) n constraints, for all 1 < A: < n : 


j)p{ip) cos^ dip = pS?.,i/p* (23) 

j)p{ip) sin((p) cos^“^((p) dip = (24) 
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We now use the maximum entropy method for deriving a ’best guess’ for the DODF that is 
consistent with this set of constraints. To this end, we maximize the entropy and use the standard 
method of Lagrangian multipliers to account for the constraints. With Lagrangian multipliers A* 
and corresponding to the constraints given by Eqs. (23) and (24), respectively, we obtain 

/ / n n \ 

p{<p) ln{p{(f)) - Xop{p) - Xip{(p) cos^{(p) - X'iP{<p) sin((^) j d(/? = 0 (25) 

\ i=l i=l ) 

C n n \ 

ln{p{ip)) + 1 - Ao - X! cos"-((/?) sin((/?) cos"-“^(v9) | 6(p = 0. (26) 

i=l i=l J 

Since this must be zero for any 6(p, the terms inside the parentheses must add to zero, which gives 


p{ip) =exp -/r-^AiCos*((/?) -^A'cos'" ^((p)sin((p) j . 

\ i=l i=l / 

The Aj and A( are found when the constraints are satisfied. The first constraint gives: 


(27) 


1 = 


^p{p)dip = e ^^exp ^ AiCos*(y?) + A'cos"' ^(y^) sin((/j)^ (28) 


The integral is called the partition function of the distribution, 

:cos*((/?) + A( cos"'“^((/?) sin((y9) 


Z(Ai, ...,An) /exp | - Y 

\ 2 = 1 


Hence 


p = ln(Z(Ai,..., A„)) 

The remaining series of constraints are satisfied by 


(29) 


(30) 


(k) 

Pii^i 


(k) 

P 21...1 


= j) p{(p) cos''{(f) dup = e '^^cos®((/9)exp ^ AjCos®((/9) + A'ces"" ^(v9)sin(99) 


1 5z(Ai,a;,...,a„,a;) d\nz 


dXi 


dXi 


(31) 


= <j>p{p) cos'{p) dp = e /cos®((/3) exp ( — ^ Ai cos®((/3) + A( cos"" ^((/?) sin(:/5) 

\ i=l 


1 5z(Ai,a;,...,a„,a;) ^ d\oz 

z dK OK 


(32) 


which gives a system of 2n implicit equations. Solving this system of equations we can find the 
2n unknown Lagrangian multipliers Aj and A( and determine the corresponding maximum-entropy 
estimate of the DODF: 


P{^) = ^ exp 


Y, cos'{p) -|- A' cos” ^{p) sin((/7) 
2=1 


(33) 
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To understand the relevance of this function we need to discuss it from an information-theoretical 
point of view. The more we know about the orientations of dislocations, the smaller is the infor¬ 
mation actually contained in the variable ip. If we take the extreme case where all dislocations are 
straight parallel edges of one sign, the information contained in the variable (p is actually zero since 
determining the orientation of any randomly chosen dislocation would not tell us anything that 
we do not already know. Maximizing the information under a set of constraints is tantamount to 
choosing a DODF which accounts for the information contained in the constraints - which reduces 
the information gained by determining the orientation of a dislocation - but no other information. 
Now, if the constraints are the known variables of a dislocation field theory, the maximum entropy 
principle allows us to construct a DODF that is consistent with the values of these variables but 
does not imply any knowledge about the orientation distribution of dislocations that is not con¬ 
tained in the given theory. Conversely, any other choice of DODF would imply that we pretend to 
know more about dislocation orientations than we can actually learn from our theory. The MIEP as 
used here provides a systematic method for dealing with the incompleteness of information that is 
implicit in any simplified, statistically averaged theory: It provides mathematical tools for handling 
information. The Shannon entropy invoked in this principle is a generic measure of information. 
Even though in certain cases a thermodynamic interpretation of the Shannon entropy is feasible 
or even illuminating, we emphasize that in the present case the information entropy associated 
with dislocation orientation distributions has little to do with the thermodynamic entropy of the 
dislocation system. We would, in fact, strongly object to any interpretation that might lead the 
reader to the conclusion that the properties of dislocation systems might be deduced from a generic 
thermodynamic principle of entropy maximization. 


In practical terms, we can use the DODE given by Eq. (33) for two purposes: Eirst, it allows us 


to provide an approximation to the exact DODE based on a limited number of alignment tensors. 
We can compare this approximation with the exact DODE that we may derive either from hdCDD 
or from systematic analysis of DDD simulations, and use this comparison to assess the performance 
of a given ODD theory. This is illustrated below in Section 4.4 Second and more importantly, we 


can use the maximum-entropy DODE, Eq. (33), to evaluate dislocation density alignment tensors 
of order higher than n. Since the parameters of the DODE depend on the alignment tensors 
of order up to n, this provides us with a method of closing the infinite hierarchy of equations, 
Eq. (16), at any desired order. We note that the resulting closure relations are optimal in the 
information-theoretical sense: Any other choice of closure relation, such as the ad-hoc relations 
used e.g. by Hochrainer et al. (2014); Hochrainer (2015), implies a different form of the DODE 


with less information, and is therefore tantamount to introducing hidden assumptions about the 
orientation distribution of dislocations that are not covered by the information contained in a nth 
order ODD theory. 

In order to use the DODE (331, we need to evaluate the Lagrangian multipliers P 21 ... 1 )) 


KiPu in terms of the alignment tensor components. The analytical solution of the system 

of equations is, except in the case n = 1 studied by [Monavari et ah] ( |2014[ ), not straightfor¬ 
ward. Instead we use a numerical approach where we tabulate the functions i)> 

Using these tabulated functions then allows us to evaluate alignment tensors of 
order higher than n from the corresponding DODE, 


= pVfc(Ai, a;, ..., A„, A'J = ..., 


(34) 


Closure at first order requires consideration of the case n = 1,A: = 2. Closure at order n = 2 
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requires to evaluate these functions for n = 2,k = 3. In both cases which are studied explicitly 
later on we find that the functions f and g can be well approximated in terms of elementary 
functions which provide us with explicit semi-analytical closure relations. 


3.5. Symmetric DODF and dislocation moment functions 

Assuming the DODF to possess an axis of symmetry w.r.t. some angle ipp can further simplify 
the closure assumptions. To understand the implications of such an assumption we note that 
anisotropy of the dislocation orientations may arise from two sources: (i) the dislocation velocity 
may be anisotropic because of some external constraint - a paradigmatic example being dislocations 
piling up against an internal boundary; (ii) anisotropy may arise from the fact that edge and screw 
dislocations may have different properties (line energy, mobility, ability to cross slip....). If (i) is the 
only source of anisotropy, then the DODF will be symmetric with respect to the angle which marks 
the orientation of the boundary. If (ii) is the only source of anisotropy, then the DODF will be 
symmetric with respect to the screw and edge orientations. Both effects may work in conjunction. 


such as in the piling up of edge dislocations against the walls of a PSB microstructure (Mughrabi 


19831 where the DODF must be symmetric with respect to the screw orientation. Only if the 
’external’ source of anisotropy (i) and the ’internal’ source of anisotropy (ii) favor different axes 
of symmetry, as in the example of a persistent slip band impinging obliquely on a grain boundary, 
the assumption of a symmetrical DODF is unwarranted. 

If the DODF is symmetrical, the alignment tensor of GND density (the dislocation density 
vector) points along the direction given by the symmetry angle ipp. The corresponding line direction 
and its perpendicular are s 


IP = p(h/|p(i)| = = [cos((pp),sin((y9p)] and = l^e = [sin((y9p),-cos((^p)]. (35) 


Transforming the coordinates of the orientation space to = ip — (p p gives: 

pi{n) p* ^p(^);)Z(^)®”'dV’ 

where Z(p) = [cos((p - (pp),sin((p - ipp)] 

The first components of the tensors are related to the moment functions: 

^p(n) _ Pii...i(^) _ ^p(.0) cos"'('(/^) = ^p{ip) cos"'((p — (pp) dip. 


(36) 

(37) 


(38) 


All other components of p'^”^ are either zero or can be calculated from the series. Dislo¬ 

cation moment functions characterize the main features of DODF. For example the first moment 
is the ratio of GND density to total density and the second moment describes the orientation of 
SSD w.r.t GND density: 


= \p^^^\/p\ (39) 

= [pfpK + 2 PS 1 K + p'SiK) Ip'- (40) 

Fig. [^depicts the different DODF constructed by the first two moment functions. 
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Figure 1: Symmetric DODFs constructed from the first two moment functions and 

plotted is 2'np{ip) as function of the angle </?,</? = 0, tt are screw orientations, ip = 7r/2, 27r/2 are edge 
orientations. From left to right: DODF with isotropically distributed SSD, GND dominated DODF 
with GND in screw orientation, DODF dominated by screw dislocation dipoles without GND, and 
DODF with both geometrically necessary screw dislocations and edge dislocation dipoles. While 
these DODFs are admissible in GDD^^\ only the first two DODFs can be represented by GDD^-^^ 


For symmetrical orientation distribution functions, the first n terms of the dislocation moment 
function series together carry the same information as the first n terms of the dislocation 

alignment tensor series In fact, both series can be converted into each other: 

^ ^ cos”"((/?) sin"((^) dip 

m times n times 


= p{r, ip) cos™(V' + ipp) sin”(V' + ipp) d<y9 

= p{r, ip){cos{'tp) cos{ipp) — sin('0) sin((/9p))™'(cos(V^) sin((/9p) + sin(V’) cos{ipp))"‘ dip 
= fpir^ip) £ cos'^~\'ip)cos'^~\ipp)sm\^p){-sm{ipp)y 


^ W j cos"" ■^(?/l)sin” ^{ipp)sm^{'ijj)cos^{ipp) 


j=o 


dip 


i=0 j=0 


^ I ( ^ i—j 


Cosm+n-i-i(^) sin*+i(^) 


X COS™- *((/?p)(-sin((^p))*sin"- ^(ipp) cos^(ipp) dip 


(41) 


where (™) denotes the number of possibilities to choose i elements out of a set of m. After removing 
those terms which vanish because of symmetry, (41) can be reformulated in tensor notation as 


k times ^2/ 

_✓'s_^ /■ ^ -N 

p(^)/pt = ®IP^ IP (42) 

(4) are F-L 

/--s 

+ IP^01P^0---01P + ... 
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Following the procedure from Section 3.4 the symmetric DODF and its moments become: 


P{^) = ^ exp ( - ^ Ai cos\ip) 


2 = 1 


^ ^ exp ^ \i cos*(V’)^ cos^('0) d'i/i, 


(43) 

(44) 


and again the unknown moments are functions of Lagrangian multipliers or known dislocation 
moments: 


= /(Ai,..., A„) = 5 (m(i), ..., m(")) = /i(p(i),..., pW). 


(45) 


3.6. Derivation of the lowest order simplified eontinuum disloeation dynamies theory - 

The lowest order simplified theory (CDD^^^) is constructed by only considering the total dis¬ 
location density p* = the dislocation density vector p(^), and the total curvature density 
g,t = q(o). The evolution equations for p*, and g* as derived from (15|-(18| then take the form 


dtp^ = V • {ve ■ p^^^) -|- vq^ 

= V • {vsp^) 

dtq^ = V • (uQ(^) - p(2) . Vu) 


(46) 

(47) 

(48) 


For this theory, the closure problem is to express the quantities and p^^'l in terms of q^,p^, 
and p^^\ Exact solutions for this problem are available in two limiting cases: (i) in the limit of 
isotropic dislocation arrangements without geometrically necessary dislocations, the second order 
alignment tensor is p^‘^'1 = ^ where is the second rank unit tensor; (ii) in the limit where 

all dislocations are geometrically necessary, we find that p^‘^'1 = \p^^'l\(l^0l^). Any suitable closure 


assumption should satisfy these two limits. Following the procedure in Section 3.5, (271 for the 
case of CDD^-^^gives the DODF as: 


pitp) = exp (—p — Ai cos(V')). 


(49) 


p and Ai are the unknown Lagrangian multipliers and are subject to solution of ( |31[ ). (42) gives 
the unknown alignment tensor p^^^ in terms of dislocations moment functions 


(2) = pt 0ip + {i- 0 ip^] , 


(50) 


and are smooth functions of Lagrangian multiplier Ai and therefore can be ap¬ 
proximated by the known moment 


M ( 2 ) ^[2 + ( M (^))2 -h (M(^))®]/4, = |p(^V/0* 


(51) 


In the derivation of the lowest order ODD (simplified ODD) in (Hochrainer et al., 20141 instead 
a linear interpolation between the two limiting cases was used: 


.(2) ~ 1 


2 L 


(/ + |pW|)Z^®Z^ + (p' 






(52) 
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Interestingly this assumption resembles the structure of (50). In the limits of = 0 (isotropic 


dislocation arrangement), or = p* (fully polarized dislocation arrangement), both expressions 
are equivalent and become exact. However, the derivation based on the MIEP provides a nonlinear 
interpolation which in the regime of small p^^^ differs significantly from the linear approximation. 
An obvious drawback of the linear interpolation is that it is non-analytic in the point p^^^ = 0 
(Sandfeld and Po 2015). A further closure approximation is needed for the curvature density 


vector. An expression which covers both limit cases is given by 

gW = -(pW)^4- 


(53) 


This expression corresponds to an equi-convex dislocation microstructure where all dislocations 
(irrespective of orientation) share the same curvature. Alternatively, the curva ture vector can be 
calculated from the approximation for p(^) using (20). Monavari et al. (2014) have studied the 
consequences of these different closure assumptions. 

3. 7. Derivation of the second order simplified continuum dislocation dynamics theory , 

While the field variables of the lowest order of the CDD equations (CDD*-^^) carry the infor¬ 
mation about total dislocation density and line curvature together with geometrically necessary 
screw and edge densities, the total screw and edge dislocation density can only be estimated. The 
theory can already represent a wide range of dislocation systems where the dislocations are 
mostly geometrically necessary, or where statistically stored dislocation are isotropically distributed 
over the possible orientations. However, the theory cannot distinguish between statistically stored 
screw and edge dislocations and must, hence, fail when dislocation arrangements become strongly 
anisotropic. In order to represent dipolar dislocation distributions it is necessary to incorporate 
the second-order alignment tensor p^^^ into the formulation . Clos ing the evolution equation (16) 
at the second order alignment tensor p^^^ gives (Hochrainer 2015): 


5ip(i) = V • {ve ®W(p(2))) 
dtp^‘^'^ = V • {vs (g) p^^^) -|- — e • p® • Vv 

dtq^ = V • (uQ(^) - p(2) • Vv) 


sym 


(54) 

(55) 

(56) 


Thus, in order to close the theory we need to represent the p^^^ and alignment tensors in 
terms of the lower-order tensors contained in the theory. Following the same procedure as before, 
we use the MIEP to reconstruct the DODF 


p{Tp) = exp p — Al cos{ip) — X 2 cos^ 


(57) 


and use tabulation to relate the Lagrange multipliers Aip to the Moment functions Fig. 

shows that the functions and are smooth non-singular functions of the Lagrangian 

multipliers. Fig. [^illustrates the shape of the DODF given by (57) for different values of the moment 
functions. 

We can then evaluate the p^^^ tensor from Eq. (42) as 


(3)/p* = {IP 0 IP^ (g) IP^ + IP^ <^IP<^ IP^ + IP^ 


I IP^ 0 IP) 
(58) 


16 





















Figure 2: Visualization of the moment functions as functions of the two Lagrangian 

multipliers Ai^2 of a second-order MIEP orientation distribution function. The red line corresponds 
to A2 = 0 which is the subspace of DODFs resulting from CDD^^^. 
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Figure 3: Left: as estimated from mO) and using the reconstructed DODF of 

CDD^ltheory. Center: Its analytical approximation: ~ MO) Vright graph: Error 

of the analytical approximation for Red lines: the admissible moment in cddO). 


In case of a pure GND microstructure this equation reduces to ^ p(i) (g) 

which is Hochrainer’s ad-hoc assumption for the tensor. We instead use the full expression 
where we approximate the MIEP estimate of as ~ MO) \/Ml^). This approximation is 
shown in comparison with the MIEP result and the approximation error in Eig. 

The second order auxiliary curvature-density tensor can be derived based on which 
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from (20) is given as the divergence of the tensoi0 


q(2)= 


with 


( 2 ) _ ^ 


qv ^ = 


2 L 


divp(3) 


J sym 


(59) 


One can also construct from the information provided by q* and following the same steps 
as for the derivation of p^^^ in the theory: 


Q(2) = 


2|q(b|2 _ 


(1 + ® + (1 - 


.(1)- 


(60) 


where Hochrainer (20151 makes again a linear assumption by setting <1> = Based on the 

MIEP we obtain instead <1> « (|q^^^|/q*)^(1 + (|q|/<?*)^)/2. 

3.8. Numerical implementation 

For the numerical implementation of our ODD models we reformulate the evolution equations 
(15) and (16) and (18) in a conservative form which is the most suitable representation of the 
underlying physics and e.g. has the advantage of numerically preserving the total number of 
dislocations in the computational domain (Ebrahimi et al.[ [2014 1: 


dtVii = V • Fi(u) + Gi(u) 


(61) 


The ‘container variable vector’ u, the vector of test functions w, the flux matrix F(u) and the 
source vector G(u) are defined as 



p* 


u = 

pin) 

, F(u) = 


qin) 


ve ■ 

vs ® 

^Q(^+1) _ £ . p{n+2) . 


(62) 



+ 


r 1 


vq 


Wp 

G(u) = 

(n — l)vQ^^'i — (n — l)e ■ p(”+^) • Vu 

, w = 

(n) 

W}, 


0 


(n) 

Wq 


For numerically solving this system of equations we use a semi-discontinuous finite element method. 
In the finite element space the weak formulation of the CDD conservation law is obtained from 
integration by parts against the vector of suitable test functions w such that 

((9tUi,Wi)o = (V • Fi(u),Wi)n -h (Gj(u),Wi) -h (Hi(u+, u“, n), ’w+)gf 7 (63) 

Ri (Fi(u), V • Wi)n -h (Gj(u),Wi) -h (Hi(u+,u“,n), w+)an + (cVu, Vw), (64) 

where (a, 6)o = We impose a periodic flux boundary condition using the Lax-Friedrich 

numerical flux H(u~*',u~, n) = |(F(u^) •n-|-F(u~) •n-|-Q;(u^ —u“)). The superscripts and 
denote the inner and outer traces of a function w.r.t. a finite element cell. The diffusive term cAu is 
introduced solely for purposes of numerical stability. The explicit third order Runge-Kutta method 


^Note that this analytically exact equation demands a numerical implementation of equations which can provide 
the gradients terms of E.g. in our Galerkin finite element framework, which is (7° continuous, this can be 

achieved by solving for as additional variables. 
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is used for time integration which benefits from both low memory usage and high stability proper¬ 
ties. This system of partial differential equations is implemented using the Differential Equations 


Analysis Library (deal.II) (Bangerth et ah, 20151. In the following benchmarking exercise where 


we also consider alternative models proposed in the literature, the same numerical implementation 
was - with small adjustments - used for all considered continuum models. This guarantees that 
deviations in the results are not due to different strategies for numerical implementation. 


4. Numerical validation and benchmark test 

For numerical validation of our kinematic closure procedure we assess the ability of hdCDD, 
and to represent the evolution of dislocation microstructure in a given velocity field, 

using discrete dislocation dynamics simulation as a reference. We also benchmark the performance 


of our equations against the continuum models proposed by Groma (1997| and Arsenlis et al. 


(20041 as well as by Acharya and Roy (20061 and Fressengeas et al. (2005). 


A good benchmark test should allow to interpret as many details of the microstructure evolution 
as possible in physical terms, and it should also allow to demonstrate where the limitations of one 
or the other model are. At the same time, it is desirable to only test a limited number of model 
aspects, in order to avoid a large numbers of parameters and interdependencies. We design the 
benchmark such that a suitable theory needs to be able to 

• differentiate between the kinematics of loops and straight lines, 

• handle the mutual conversion of SSDs and GNDs; we consider this as an essential pre-requisite 
for any averaged theory that operates on scales above the resolution of individual dislocation 
lines, 

• differentiate between SSDs of screw or edge orientation. 

An ideal benchmark situation in this respect is the motion of mobile dislocations in a persistent slip 


band (PSB) microstructure Fig. 4aI, for which abundant experimental evidence is available (see 
e.g. [Mughrabi (1983|) including in-situ observations of dislocation motion (see e.g. Lepinoux and 


Kubin (19851), such that the simulation results can be directly related to experimental information. 


4 . 1 . Formulation of the benchmark problem 

As benchmark system we consider a PSB structure consisting of a periodic arrangement of 


obstacle-rich walls separated by obstacle-free ’channels’ (see Fig. 4b I. The ’walls’ mainly consist of 


immobile edge dipoles which hinder the motion of dislocations but do not themselves participate in 
the plastic deformation process, while in the ’channels’ dislocations can move freely. Note that we 
do not consider the formation of such a persistent structure during cyclic loading (for an overview 
of theoretical models for PSB patterning the reader may consult e.g. Kubin et al. ( 2002| )). Here, 
we only pursue the modest aim of modeling how dislocations move in such a strongly heterogeneous 
microstructure by investigating the motion of dislocations over half a loading cycle, in a situation 
where the obstacle microstructure has already reached a quasi-stationary state - which may persist 
over thousands of cycles. We describe dislocation-obstacle interactions in terms of a yield stress 
which depends on the density of immobile edge dislocation dipoles in the walls, = aqib 
where p, is the shear modulus, b the Burgers vector length, and a ~ 0.2 a non-dimensional coefficient 
characteristic of dipole hardening. The yield stress reduces the resolved shear stress such that 
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the velocity is given by f — r^) where B is the dislocation drag coefficient. In the 

following we assume parameters typical of PSB microstructures in Cu, viz, B = 1 x IQ-^Pa s, 
b = 0.256 nm, /r = 4 x 10^ MPa, a = 0.2, and = 6.2^^ m“^, corresponding to a typical stress 
level for wall deformation of ~ 60 MPa. The dipole density drops from in the walls to zero 
in the channels, and we model the smooth transition in between by a sigmoidal function. The 


concomitant velocity profile is shown in Fig. 4c and assumed time independent. The assumption 
of a stationary velocity field can be justified by looking at experimental data: The density of 
immobile dipoles in a fully developed PSB structure (10^^ — 10^® m“^) is about two orders of 
magnitude higher than the density of the moving dislocations (10^^ — 10^^ 


m 


-2 


see 


Essmann and 


Mughrabi (19791; Mughrabi (19831. This has two consequences: (i) the dipole density is a slow 
variable which cannot change much over a half cycle, hence we can assume the obstacle field and 
the concomitant velocity field as stationary, (ii) the flow stress is controlled by interactions between 
mobile dislocations and immobile dipoles, rather than among mobile dislocations. Hence we can 
assume the dislocation velocity to be independent of the densities of the moving dislocations. 

In our DDD reference model we consider a simulation volume of size La, x Ly x = 1.4 x 15 x 
20 pm^ with periodic boundary conditions in all three spatial directions. The simulation volume 
contains 900 mobile dislocation loops which are initially assumed circular with radius rg = 0.4/im 
and which are located at random positions. This configuration corresponds to an initial mobile 
dislocation density p* = 5.3 x 10^^ m“^ = pg. The large values of Ly and Lz are chosen in order to 
have sufficiently many loops for obtaining meaningful statistical averages. 


J^.2. Brief formulation of the benehmark theories 

In the following we benchmark the performance of CDD and of other models proposed in the 


literature. We consider the model of Groma (1997| for straight parallel dislocations, the edge-screw 


dynamics of 

Arsenlis et al. 

(2004 

1, and two versions of the semi-phenomenological model of Acharya 

and co-workers (Acharya and Roy 

2006 

Pressengeas et al. 

20051. 


The problem we consider is statistically homogeneous in the y and z directions and the corre¬ 
sponding simulations can be envisaged as simulations of ”1.5D” systems: the computational domain 
is one-dimensional, dislocation densities depend only on the x coordinate, and any geometrically 
necessary dislocations must possess edge orientation. However, at the same time the evolution of 
dislocation densities is influenced by two-dimensional dislocation motions. With the exception of 
the model of Groma, which considers straight parallel edges only, all considered models contain 
information about ’perpendicular’ motion of dislocations in y direction. Some of the models con¬ 
tain such information in explicit form through variables such as the curvature density in GDD^^^ or 
gdd( 2) , or the equations describing screw dislocation motion in the model of Arsenlis et ah. Oth¬ 
erwise such information is implicitly contained in parameters describing dislocation multiplication, 
which is a process that necessarily requires two-dimensional dislocation motion. 

4-2.1. The referenee DDD model 

DDD simulations serve as reference providing detailed insight into the discrete microstructure 
evolution. The used DDD model is based on the so-called Parametric Dislocation Dynamics method 
(PDD), cf. (Ghoniem et al. 20001 as implemented in Po and Ghoniem| ( 2014| . There, cubic 
splines are used in order to have non-vanishing curvature along each dislocation segment. Since 
we are here only interested in the kinematic consistency of different models, we evolve the discrete 
dislocation lines in the same analytically given velocity field as used in the GDD simulations. 
The evolving discrete dislocation microstructure can be directly compared with the respective 
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Figure 4: (a) Persistent slip band structure: Low density channels are separated by high density 
walls which mostly consist of edge dislocations in dipolar formation. Note that dislocations with 
the same orientation have the same curvature, (b) Sketch of an idealized PSB: dislocations of edge 
character get trap and bow-out from walls and screw segments can easily glide through channels 
. (c) The kinematic setup of idealized PSB: Prescribed velocity is maximum in the channels 

and drops to 15% in the middle of the wall (darker shaded area); the initial mobile dislocation 
population consists of uniformly distributed glide loops; the dashed line indicates the periodic 
simulation box. (d) Initial dislocations segments are distributed uniformly in all direction which 
means total dislocation density is consisting of equal amount of edge and screw components. 


continuous fields through the Discrete-To-Continuum (D2C) strategy (Sandfeld and Po, 20151: for 
extracting continuous field data from DDD we start from the discrete line description and convolute 
them with a non-local smoothing kernel function to obtain continuous and differentiable continuum 
data that can easily be compared with data obtained from hdCDD/CDD. Further details on the 


D2C strategy can be found in (Sandfeld and Po 2015). 
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4-2.2. The model of Groma: straight positive and negative edge disloeations 

Groma and co-workers (Groma 1997 Groma et al.[ |2003 | developed a kinematic theory for 
systems of straight parallel edge dislocations moving in ±x direction. The microstructure is divided 
into positive and negative p~ edge densities (Fig. which evolve according to 


dtp^ = -dx (^p~^vj and dtp = dx {p v) , 


Total density and edge component of GND vector follow from (65) as 

p^ = p'^ + p~ and 


= P'^ 


P 


(65) 


( 66 ) 


To be consistent with our initial conditions, we choose the initial conditions for the Groma model 
as p* = po,p^ 2 ^ = 0 corresponding to p'^{x) = p~{x) = po/2. 

Note that in the Groma model all dislocations are assumed to be straight, and therefore there 
can be no dislocation multiplication (the total dislocation line length is constant). We note that 
various authors have introduced phenomenological generalizations of the model to account for dislo¬ 
cation multiplication - which together with an appropriate choice of the concomitant fit parameters 
would make the model perform much better in the following benchmark exercise. Here, however, 
we are interested in the model precisely because of its very simple and transparent form which 
does not contain any ad hoc assumptions regarding dislocation kinematics. 


4-2.3. The model of Arsenlis: serew-edge representation of dislocation kinematics 


Arsenlis et al. (2004) explicitly acknowledge the distinction between dislocation kinematics and 


dynamics. In order to make the kinematic problem computationally more tractable, they restrict 
the dislocation orientation degrees of freedom to a small number of line directions by taking into 
account only the four orientations p®"*", p®“, p®"*", p®“, where superscripts e and s denote either pure 
edge or pure screw densities, respectively, and the sign (-I-) or (—) indicates the polarity of the 
dislocation density (Fig. [^. In the 1.5D case and assuming that all dislocations of all orientations 
have the same mobility and thus the same velocity in a given point, the evolution equations reduce 
to 

dtp^+ = -dx ((p"+ + fp^)v) + q^v 

dtp^- = dx {{p^- + (1 - r)p» + q^v 

dtp" = (g"+ + q^~)v 


(67) 


P = 


pS+ ^ 


and the function (for explicit expression see Arsenlis et al. (2004)) is 


needed to ensure that dislocation densities do not become negative in the course of evolution. 
In our presentation of these equations we have emphasized the analogy with GDD theories by 
introducing the notations = p^~^/l^~^ ,q^~ = p^~/l^~, and = p®/P where l^~^, and T are 
average dislocation segment lengths . These quantities play in the Arsenlis formulation very much 
the same role in controlling dislocation multiplication as the curvature density in GDD theories. 


Accordingly, Arsenlis et al. (2004) complement the equations (67) by equations of evolution for 
and T. These equations are the counterpart of the equation for curvature density evolution in 
GDD theories. We do not give them explicitly since in the absence of dislocation sources - which 
is the case in the present benchmark problem where the number of loops is fixed - the equations 
given by Arsenlis et al. (2004) imply that q^~^,q^~ and q^ are time independent such that we can 


replace them by their initial values. 
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Figure 5: Schematic depiction of the DODF implied by the Groma, Arsenlis and Acharya models. 
Note that the peaks symbolize Dirac delta functions. Left: Groma ± uses edge dislocation densities 
of both signs, center: Arsenlis describes the dislocation microstructure with four densities of ± 
edge and screw orientations, right: In Acharya’s model the dislocations are divided into a GND 
density with known direction and SSDs whose orientation distribution is unspecified. 


To be consistent with the initial values of the dislocation density and dislocation density vector 
in our DDD simulations we set the initial dislocation densities as /2 = po/4. Furthermore, 

to correctly reproduce the initial dislocation line length increase we set = q^ /2 = p*/(4ro). 
For purposes of comparison with the other theories we consider the total dislocation density p* = 
p*^+ + p®“ + p®, the edge component of the GND dislocation density vector, = p®"*" — p®“, and 
the mobile SSD density p®®° = P* — 


4.2.4- Models of Acharya and co-workers: Combining the dislocation density tensor with phe¬ 
nomenological models for statistical dislocation density evolution 

(20061 starts out from the classical kinematic 


The crystal plasticity model of Acharya and Roy 


equation for GND transport formulated by Mura (19631, dta = —V x [a x v) where v is the 
velocity vector of the geometrically necessary dislocations. In order to account for the contribution 
Lp of statistically stored dislocations to the plastic distortion rate, this is generalized to 


dta = —V X (a X u + Lp), 


( 68 ) 


where Lp is the plastic distortion rate due to the motion of SSDs. Acharya and co-workers do not 
make attempts to derive the temporal evolution of Lp from the underlying motion of dislocations, 
but rather consider this as a problem which should be solved in the spirit of classical continuum 
mechanics by making constitutive assumptions on a phenomenological basis. 

We consider two cases: (i) the constitutive assumptions made by Acharya and Roy (20061 
(henceforth referred to as ’Acharya’) suggest to set 


Acharya and Roy (2006 


Lp = Pp^bv. (69) 

where P = {b®n)/b is the slip system projection tensor and p™ is implicitly assumed to account 
for the statistically stored, mobile dislocations only. Further assumptions in the same paper are 
tantamount to using the same velocity modulus v for geometrically necessary and statistically 
stored dislocations (very reasonable given that any individual dislocation has no means of decid¬ 
ing whether it is ’geometrically necessary’ or ’statistically stored’) and considering the density of 
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statistically stored mobile dislocations as constant. 


Fressengeas, Varadhan and Beaudoin (Fressengeas et al. (20051, henceforth referred to as FVB) 
instead proposed to derive the evolution of the mobile SSD density from the phenomenological 
Kubin-Estrin model (see e.g. [Kubin et al.| ( |200^ ) which yields: 


dtP^ = [Ci/b^ - C2P“ - 


(70) 


where Ci accounts for mobile dislocation multiplication, C 2 accounts for mobile dislocation an¬ 
nihilation and/or dipole formation, and the term with C 3 accounts for immobilization by forest 
obstacles. 7 is the total plastic shear strain rate in the slip system and includes both SSD and 
GND contributions. Since forest obstacles are absent in our benchmark example which considers 
dislocations on a single slip system, we only consider the first two terms. (Note that considering a 
dipole formation term is not inconsistent with our assumption of a quasi-stationary arrangement 
of dipole-like obstacles, since in actual PSB microstructures the dipole densities exceed the mobile 
dislocation densities by about two orders of magnitude and change little over one loading cycle). 
The parameters Ci and C 2 are fitted bona fide to reach, at the end of the simulation, the same 
average plastic strain and mobile dislocation density as found in our benchmark hdCDD theory. 

In both cases (i) and (ii) a corresponds to via cx = p^^^ ® b and the total density p* is 
equal to p™ -|- |p^^^|. Accordingly we set the initial values p™ = po,Q: = 0. For comparison with 
the other theories we consider the total dislocation density p*, the edge component of the GND 
dislocation density vector, = a 2 ilb, and also the mobile SSD density p™. 


4-3. Results from DDD 

Reference data are obtained from DDD simulation with an initially random distribution of 
loops. Fig. shows typical dislocation microstructures at three different time steps. Note that the 
loop configuration is projected along the 2 ; axis into the xp-plane: The different loops expand on 
different slip planes. A number of important features can be observed: dislocation loops in the 
channel expand until they reach the wall. There, line segments change their orientation such that 
they align with the contour lines of the velocity field. Parts of dislocation loops in the channel 
regions further expand until they reach near-screw orientation and become almost straight. Since 
the wall regions are not completely impenetrable, dislocations eventually cross the wall and leave 
on the other side to bow out in the adjacent channel. A significant number of dislocations is getting 
’trapped’ inside the walls and form near-parallel edge dislocation bundles of zero net Burgers vector. 
In the transition region between the wall and the channel we observe that the radius of curvature, 
with which lines are bent, is smaller than the initial loop radius. 

4.4- Results from hdCDD and vs. DDD: analysis of disloeation density in the eonfigu- 

ration spaee 

For a detailed comparison of the DDD microstructure with GDD results we convert the discrete 
lines into a continuous density in the higher-dimensional hdGDD configuration space, here spanned 
by the variables x and <p. This offers the possibility to directly compare with hdGDD, and with the 
reconstructed angular densities derived from GDD*-^^ and GDD^^^. Since the benchmark system 
is statistically homogeneous in y and 2 ; directions, the angular densities depend on x only. We 
obtain continuous data by first applying the D2C conversion to the three-dimensional DDD data. 
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Figure 6: Evolution of an initially random distribution of discrete dislocation loops in the ’PSB’ 
velocity field shown in Fig. . From left to right: initial values and dislocation microstructure at 
two subsequent time steps, corresponding to strains of 0.33% and 1%. 


followed by averaging over the y and directions. Higher-dimensional density distributions which 
result from this process at strains of 0.33% and 1% are shown in the top row of Fig. 

In principle it is possible to use the spatially averaged initial DDD microstructure directly as 
initial condition for the continuum simulations (compare Sandfeld and Po 2015[). Instead, we 


consider as the initial state the statistically homogeneous and isotropic densities which correspond 
to an average over the initial conditions in an infinite ensemble of DDD simulations, or to the 
spatial average obtained for a system of infinite extension in y and z directions. This leads to the 
following initial conditions: The hdCDD dislocation density and curvature density functions are 
isotropic and spatially homogeneous, p{r,ip) = pQ/{2 t:), q{r, ip) = /9o/(2'?rro). The corresponding 
initial values for the CDD^^^ and CDD^^^ variables are p* = poj = Po/^O) = 0, p*-^^ = 

For hdCDD, the evolved density function p(r, ip) can be directly compared to the DDD data. For 
CDD^^^ and CDD^^^ we evaluate the DODF from the field variables in each spatial point using the 
maximum entropy principle as discussed in Section 3.7 and then multiply with p* to obtain the 
matching density function p(r,p). 

At first dislocation loops expand uniformly in the channels (x ~ ±(300 ... 750) nm), where in 
the hdCDD configuration space positive edge orientations (p = O.Stt) move towards the right and 
negative edge orientations (p = l.Svr) move towards the left. Eventually more and more density 
from the channel reaches the wall where positive and negative edge components of dislocations 
accumulate on the left and right sides of the low velocity region, respectively. These accumulations 
show in Fig. [^as spots of high hdCDD density. We observe that both hdCDD and CDD^^^ are in 
very good agreement with the reference DDD simulations in terms of the morphology and location 
of the density accumulations. CDD^^^ also shows accumulation of GND in edge orientation at the 


25 













2.()7r 


7 = 0.33% 


7 = 1% 



2.()7r 






. 4 . - 




f : 



.. . . J . 





-600-400-200 0 200 400 600 

X [iiml 


-600-400-200 0 200 400 600 

X fninl 


J L. 


-I 


9 13 18 22 26 31 35 

[ 10^2 771 “^] 


0 9 18 26 35 44 52 61 70 

p^ [ 10 ^^ 777 “^] 


Figure 7: hdCDD density p{x,(p) in the configuration space: the vertical axis denotes the line 
orientation, the horizontal axis is the x-axis; top row: p(x,ip) as reconstructed from DDD, second 
row: p{x, ip) as obtained by solution of the hdCDD equations, third and fourth rows: p{x, p) as 
reconstructed from the solutions of the CDD^^^and equations. 
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walls but fails to capture the concomitant transition to screw orientations in the channels. 

At a strain of 0.0033 the channels are already partially depleted of edge dislocations which 
accumulate at the left and the right edge of the wall whereas the dislocation density at the center 
of the wall (around x = 0 ) is still low because the small velocity makes it difficult for dislocations 
to penetrate the wall. This is well captured by all models. Dislocation loop segments approaching 
the wall assume more and more straight edge-like orientations, whereas segments moving in the 
channels are preferentially screw oriented. hdCDD is able to represent this feature exactly, and 
CDD*-^^ shows qualitative agreement, while does not indicate any enhancement of near¬ 

screw orientations in the channels. 

When the system further evolves (strain 0.01), edges align inside the wall into dipole bundles. 
hdCDD correctly represents this through the two high density accumulations which now overlap 
in the spatial projection. In accordance with Fig. |^a) and the DDD plot in Fig. |^the curvature 
of these dislocations is almost zero. At the same time in the channel the density consists almost 
exclusively of threading dislocations in positive and negative near-screw orientations. Both the 
DDD data and the hdCDD simulations show this feature which is also correctly represented in the 
CDD*-^^ simulations. CDD^^^ by contrast, makes the erroneous prediction that the channels are 
depleted of edge and screw dislocations. 

We conclude that CDD^^^ is not able to fully predict the microstructure evolution. The seg¬ 
regation of edge dislocations in the walls and screw dislocations in the channels is not captured. 
cdd( 2 ) on the other hand, is able to represent most important characteristics of the PSB-like 
system. Only small details show deviations from DDD. Finally, hdCDD is able to reproduce all 
relevant mechanisms and the time evolution of dislocation microstructure in very good agreement 
with DDD. This confirms the observation made in previous works that hdCDD is kinematically 
exact, i.e., its results coincide with averages over an infinite ensemble of matching DDD simulations 
carried out assuming the same velocity field. 

DDD, CDD^^^and hdCDD predictions are in a good qualitative agreement with observations 
of dislocation microstructure and dislocation motion in PSB in copper using transmission electron 


microscopy (TEM) (Fig. 4aI. We quote Mughrabi et al. (19791 who stated that “(1) edge dislo¬ 
cations bow out of the walls and become entangled in the neighboring walls; ( 2 ) edge dislocations 
that are newly formed in the wake of moving screw dislocations are incorporated into the wall”. We 
note that Mughrabi ( 1987| also reports the observation that dislocations with the same orientation, 
located at equivalent positions in the wall-channel structure of a PSB, have the same curvature - 
which is a main assumption for curvature-density fields in the hdCDD and CDD theories. Also the 
observation that this curvature is high near the wall edges but lower in the channels is in line with 
the results of our simulations, see Fig. 

^.5. Comparison of results from CDD^'^^ with other theories: disloeation density/strain profiles and 
average behaviour 

To compare CDD With other theories we need to specify a minimum set of variables which 
can be found in all of the theories considered. In what follows we look at the spatio-temporal 
evolution of the total dislocation density p*, the edge component of the GND vector and the 
plastic slip 7 . We also include a derivative quantity, namely the mobile SSD density which can be 
evaluated as = p* — These quantities are shown in Fig. [s for the same two snapshots in 

time as in the previous section. For hdCDD the integrated densities p*, p^'^ are obtained from the 


higher-dimensional p(r, p) using (10)-(13), CDD^^^ already contains these field variables in explicit 
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form, whereas from the other theories and can be deduced as discussed in the respective 
sections. 

We observe in Fig. |^e) that the initial increase in GND density at the wall boundaries is 
predicted by all models. Models without dislocation multiplication (Groma) or which associate 
dislocation density increase only with GNDs (Acharya), however, predict a slightly lower GND 
density concomitant with a lower strain. Gomparing the edge GND density p^^ with the total 
density (Fig. [^a)) one observes that at the boundaries of the wall all dislocations have edge 
character which is also predicted by all models. Most of the plastic strain is generated inside the 
channel where the average dislocation velocity is higher. 

However, as the systems evolve further, only GDD*-^^ and the Arsenlis model show good qual¬ 
itative agreement with the hdGDD or DDD results (Fig. |^b)), while models without dislocation 
multiplication (Groma, Acharya) strongly underestimate the strain in the channel and the disloca¬ 
tion density in the wall. An intermediate position is provided by the FVB model with dislocation 
multiplication. This reproduces accurately the strain and GND density profiles but strongly under¬ 
estimates the density of SSDs in the wall - where it predicts a minimum rather than a maximum 
of the SSD density. This feature is important because it makes the model a poor candidate for 
explaining wall formation - which requires SSDs to form narrow dipole configurations inside the 
walls rather than in the channels. According to the FVB model with dislocation multiplication, the 
opposite is true since the mobile SSD density (the density of those dislocations which can mutually 
trap into dipoles) is highest in the center of the channel and lowest in the wall. 

We may also look at averages. All models which account for dislocation multiplication can more 
or less correctly reproduce the increase in average plastic strain, and the concomitant increase in 
the spatially averaged total dislocation density. In the case of GDD^^^ and the Arsenlis model this 
is achieved without parameter fitting, in case of the FVB model we have fitted the multiplication 
and dipole formation parameters to obtain the correct values of end strain and mobile dislocation 
density. We observe two distinct stages in the evolution of p* and 7 : in the first stage (7 < 0.33%) 
we find free expansion of loops which leads to a linear increase of the average dislocation density 
(/ 9 *) (which is proportional to the circumference of all loops) and a quadratic increase of the average 
plastic strain ( 7 ) (which is proportional to the area of all loops). In this stage most theories can 
predict the evolution of dislocation density more or less correctly. In a second stage deformation 
is mainly by screw dislocations moving in the channels while dislocation multiplication is mainly 
by the concomitant drawing out of edges along the wall. In this stage, both the strain and the 
dislocation density increase with time in an approximately linear manner. 

An instructive feature shows up when we compare the amount of dislocation multiplication 
with that in a spatially homogeneous reference system where dislocations move at an effective 
velocity Ueg = {v{x)) which is the spatial average of our homogeneous velocity field. Surpris¬ 
ingly, both hdGDD and to a lesser extent GDD*-^^predict that the heterogeneous system shows a 
larger dislocation density increase than the homogeneous reference system with uniform velocity 
(Fig. I^a)), whereas the total plastic slip generated by the movement of dislocations is less than for 
the homogeneous reference system (Fig. j^b). This is due to the following features: (i) dislocations 
accumulate in the wall where they move more slowly than on average, leading to a reduced total 
strain, (ii) dislocation segments which have crossed a wall can expand in the adjacent channel and 
effectively act as additional sources, leading to enhanced multiplication. The Arsenlis model which 
does not account for spatial transport of dislocation curvature predicts the dislocation increase in 
the heterogeneous system to be identical to the uniform reference system. This indicates that it can 
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Figure 8 : Evolution of total dislocation density p*, mobile SSD density edge GND density 

p^\ and plastic slip 7 in the slip plane. Using the (anti)symmetry of the system hdCDD and 
CDD^^^are plotted only on the left side of the domain and the other models on the right side. 
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Figure 9: Normalized total dislocation density / 0 *and total plastic slip 7 over time. 


not fully account for the effects of dislocations bowing out from a wall. The FVB model produces 
a correct description of the strain and dislocation density increase - two variables to which the two 
parameters of the model have been fitted. The Groma and Acharya models, which cannot account 
for dislocation multiplication where this leads to accumulation of statistically stored dislocations, 
by contrast, completely underestimate /?*and 7 . 


5. Conclusion 


We have applied the Maximum Information Entropy Principle (MIEP) to reconstruct the Dis¬ 
location Orientation Distribution Eunction (DODE) from a finite set of the associated alignment 
tensors and have used this reconstruction for kinematic closure of ODD evolution equations. This 
method can in principle be used for closure of CDD theories using alignment tensors of any order, 
though we have explicitly considered only the lowest two orders. 

Eor testing the performance of the resulting evolution equations, we have proposed a generic ap¬ 
proach which compares averaged DODE obtained from DDD simulations with those reconstructed 
from CDD models. We have carried out such a comparison for a benchmark example, namely 
the motion of dislocations in an idealized PSB microstructure. The comparison demonstrated the 
capabilities of boh hdCDD and CDD^^^ to account for all essential features of dislocation motion in 
such microstructures. At the same time CDD^^^ - which performs well in other situations, see e.g. 
Hochrainer et al. (2014| - was shown to be incapable of dealing with situations where edge-screw 
asymmetry is pronounced. This demonstrates that the formulation of CDD evolution equations 
requires careful consideration of the structural complexity of the problem under consideration. 

We have then benchmarked the performance of CDD^^^ against several models published in 
the literature. Out of the investigated models only CDD^^^ and the Arsenlis’ model - which is 
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specifically designed to capture screw-edge systems - were able to predict all essential characteristics 
of dislocation motion in the strongly anisotropic PBS microstructure. 

Because of its kinematic exactness, hdCDD - even though computationally inefficient as a 
simulation tool - provides a useful reference for assessing kinematic consistency of simpler theories. 
Obtaining a similar reference for assessing dynamic consistency is a more difficult task: In this case, 
we need to perform ensemble simulations of discrete dislocation systems and then use ensemble 
averaging to determine the corresponding averaged density measures and their time evolution. 
Performing an assessment of dynamic consistency of the proposed approaches, for instance for 
interacting loops under an external shear stress in the geometry outlined above, as well as for other 
benchmark problems, remains an essential task for future investigations. 
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Appendix A. Comparison of computational cost of different density-based theories 


One of the possible advantages of coarse-grained continuum dislocation theories over discrete 
dislocation dynamics is their lower computational cost. The number of degrees of freedom in a 
continuum dislocation model equals the number of state variables times the number of nodal points 
used in the representation of the density fields. Hence the number of state variables in each spatial 
point is a good criterion for comparing the computational cost of different continuum models. In 
our 1.5D benchmark test, the homogeneity of the problem in the direction of screw dislocation glide 
implies that the geometrically necessary dislocations must always have edge orientation which im¬ 
plies that = 0, p^i 2 — 0 = P^~ which reduces the number of state variables. In a fully 

3D formulation, p^^^and p^^^^remain 2D vectors and symmetric tensors in the local slip system co¬ 
ordinates and therefore have 2 and 3 components, respectively. In the higher-dimensional hdCDD 
theory, the variables p^p, have to be partitioned in the orientation space in addition to the spatial 
domain which makes hdCDD infeasible for real 3D systems. In the present benchmark problem, we 
use for hdCDD a partition of orientation space into Up = 60 regular Galerkin elements. Comparing 
the total number of DOF and state variables for each model given in Table (A.l) shows that all 
the models except hdCDD have relatively low computational cost. 



hdCDD 

cdd ( 2 ) 

cdd(i) 

Arsenlis 

Groma 

Acharya 

FVM 

variables 

Pip, qip, 7 



pe±^ps±^^ 

p-\- 

P ,1 



1.5D/2.5D 

-|- 1 

5 

4 

4 

3 

2 

3 

3D 

2np, -|- 1 

7 

5 

5 

- 

3 

4 

DOF 

~ 48000 

~ 2000 

~ 1600 

~ 1600 

~ 1200 

~ 800 

~ 1200 


Table A.l: Number of state variables per slip system at each spatial point and total degrees of 
freedom in the benchmark problem for different models. The domain is discretized with 200 second 
order elements. 


If we now compare with our DDD reference model we may calculate the number of DOF 
for the reference model as follows: In the reference model we have 900 loops, and in order to 
achieve an angular resolution comparable to hdCDD each loop needs to be discretized into 60 
segments (corresponding to an initial nodal spacing of about lOnm). This corresponds to a total of 
54000 nodal DOF. However, we emphasize that this is not a fair comparison since the continuum 
calculations use the statistical homogeneity of the benchmark problem in y and directions to 
reduce dimensionality which obviously results in a tremendous reduction of degrees of freedom - 
for a generic three-dimensional problem of the same size the number of DOF would need to be 
about three orders of magnitude higher. On the other hand, we also emphasize that the main 
differences in computational cost between DDD and coarse-grained continuum theories arise from 
the possibility, in a density-based theory, to represent segment-segment interaction stresses as local 
functionals of the dislocation densities, instead of calculating them by summation over all segment 
pairs, whereas the calculation of long range stresses can be carried out on a much coarser mesh 
than is needed for solving the transport equations. We will discuss these issues in a future work 
which presents a three-dimensional dynamical implementation that includes interaction stresses. 
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Appendix B. Dislocation annihilation in continuum dislocation dynamics 


The problem of annihilation of dislocations only emerges in coarse-grained continuum theories 
that allow for the coexistence of dislocations with different orientation within the same volume 
element. CDT theories such as Arsenlis et ah (2004| formulate this problem in a conceptual 


framework that focuses on encounters of straight lines which annihilate once they meet within 
a reaction cross-section (annihilation distance) leading to bi-molecular annihilation terms. But 
dislocations are not particles. Dislocations annihilate when two loops merge - a process which in 
itself consumes very little dislocation line length - and the subsequent motion then may lead to a 
reduction rather than an increase in line length, but at other locations and in different orientations 
from those of the annihilating segments. 

Experimental observation indicates that the annihilation of segments with near-screw orienta¬ 
tion may proceed via a cross slip mechanism and that the minimum distance ya required for this 
process to occur is significantly (by about two orders of magnitude) larger than for annihilation 
of segments of other orientations. Therefore we first consider situations where the annihilation 
process is initiated by cross slip and subsequent annihilation of two near-screw segments which are 
oriented at a small angle (fa. £ from the screw orientations (^a = 0 and ^pa. = tt. 


Fig. B.IO (left) depicts a situation some time after near-screw segments of two loops moving 


on parallel slip planes with normal vector n have merged by cross slip. As the loops expand, the 
intersection point A - which corresponds to a segment in the cross slip plane and separates segments 
of direction l((p) and in the original slip planes - moves in the Burgers vector direction to the 
new position A. During this movement of the intersection point the segment AB with orientation 


ip as well as a matching segment with orientation p' = tt — p (Fig. B.IO center) rotate into the 


respective screw orientations and annihilate. Assuming that the loops have the same curvature and 
that the dislocation velocity is independent of the segment orientation p, the annihilation process 
is symmetric with respect to the direction of the Burgers vector. 

The total annihilation rate is determined by the probability of finding, for a given segment 
of orientation p, an intersection point with a matching segment of orientation tt — p, and by the 


velocity of this intersection point parallel to the segment direction l{p). As depicted in Fig. B.IO 
(centre), if the segments of both loops are moving with velocity v, then the velocity component of 
the intersection point parallel to l{p) is 


Va. = vcot{p) 


(B.l) 


To evaluate the probability of finding for a segment of orientation p an intersection point with a 
matching segment we proceed as follows: (i) Matching segments must have orientations between 
TT — p — A(/j and tt — p + /S.p. The area density of intersection points of such segments with a 
perpendicular plane (normal vector 1{tt — p)) is given by 

^TT — 

/3a(7r -p)= p(<^')d(/7'. (B.2) 

Jir—ip—Acp 

(ii) We now consider a strip of width 2ya around the slip plane of the first segment, but still 
contained in the plane with normal vector 1{tt — p). The average distance d! between intersection 
points, evaluated along the direction of this strip, follows from the condition 2d'yaPa{'^ ~ v) — 1- 
Their spacing along the direction l{p) is then (see Fig. 

d{TT — p,p) = . = . , d' = --. (B.3) 


B.IO 


sin(, 0 ) sin( 2 (/ 5 ) 


right) 
d' = 


1 
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Figure B. 10: Top view of a cross slip annihilation process. Left: Cross slip initiates the annihilation 
of two merging dislocation loops. Gray and black lines depict the evolution of the merging loops 
during a small time step. The dashed lines shows the annihilated parts of the loops. Center: The 
intersection point A moves to a new position A'. Its relative motion to the segment (p decreases 
the segment length by the length AB. Right: segments with orientation p' have expected spacing 
d perpendicular to their line direction. If evaluated along the line direction their spacing is 

d' — 

“ “ sin(/3) • 


We now assume that the configurations of the two merging loops are statistically independent. This 
assumption implies that the loop centers are located at statistically independent points and that 
each loop expands, away from the intersection point, in a manner that is not strongly influenced 
by the presence of the second loop. In that case the volume density of such intersection points 
can be evaluated as the line length of orientation p per unit volume, divided by the mean spacing 
d{p, TT — p) of intersection points along a segment of this orientation. Accordingly, we obtain 


Ant(¥^,vr - p) 


d' 


2y^sm.{2p) 


p{'^) / 

J TT—Lp — Acf) 


(B.4) 


It follows that the temporal change of p{p) due to the motion of such intersection points (’annihi¬ 
lation rate’) becomes: 
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P{P) 

rn—ip+Acf) 

pip)dp 
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7r—(p—A<p 



pn—ip-\-A(l) 

4 

P{P) 

p{p')dp' 


J 

7r—(p—A<p 


cos((^) 

. , . s\n{2p)y^v 
sm((/?) 

cos^ {p)ya.v. 


(B.5) 

(B.6) 


We consider two limiting cases: First, if ah dislocations are screw oriented then p{p) = pj^5{p) + 
P-5 {'k — p). The annihilation rate follows as 

^iP+lanil = ^i/’-lanil = -^P+P-V^V (B.7) 

which is the result expected by kinetic theory for 2D particles moving with velocity v in opposite 
directions and annihilating if they pass within a reaction cross-section 2ya- The second limiting 
case, which is the one relevant to our benchmark example, is that the dislocations have a smooth 
angular distribution which can be approximated as constant over the small angle interval 2/S.p. 
Then, 

9tp{p)\^ni\ = -8App{p)p{Tr - p) cos^{p)y^v. (B.8) 
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We see that the screw-triggered annihilation process removes segments of all orientations with a 
rate that is largest for screw orientations and goes to zero for the pure edge orientations (/? = 7r/2 
and ip = 3'7r/2. The concomitant temporal change of the zeroth-order alignment tensor (total 
dislocation density) can be calculated by integration: 


^‘^lanil ^ I anil 


(B.9) 


The DODF of PSB structures is symmetric around the edge orientations = f and ip = As a 


consequence, p{ip) = p{tt — ip) . Using this symmetry property and the DODF given by Eq. (49), 
the rate of reduction in total dislocation density in CDD^^^can be evaluated as 


V =ipl 


t^2 8 A(/? 


anil 




exp(—2Ai sin((/?)) cos^((/?)d(/9 




(B.IO) 


where Z is the partition function of the distribution. The first order alignment tensor is not directly 
affected by the annihilation process. For completeness, we also consider the case (not relevant to 
our benchmark example) where the DODF is symmetric around the screw orientations ip = 0 and 
ip = TT such that p{Tr — ip) = p^mp). In that case an analogous calculation gives for CDD^^^ 


dtp' 


anil 




cos^(y?)d(/9 


ViiV 


(B.ll) 


For 

cases 


a completely isotropic dislocation arrangement, Ai = 0 and Z = 27r, and we obtain in both 


dtp' ={p') 


t\2 2A(^ 


anil 


TT 




(B.12) 


In the presence of excess dislocations of one sign, the annihilation rate is a function of the first 
moment function M^) = \pd-l\/p^ which can be understood as the GND fraction of the total 


dislocation density. This dependency is illustrated in Figure Fig. B.ll which shows the annihilation 
rate, normalized by the value at M^-) = 0, as a function of the GND fraction M^). 

We can see that the annihilation rate decreases monotonically with increasing GND fraction 
and goes to zero if all dislocations are GND. However, this decrease is not well described by 
the parabolic dependency expected according to kinetic theory for a system of straight parallel 
dislocations: In edge-dominated microstructures such as PSBs, where GNDs are of edge orientation 
but annihilation is triggered by cross slip of screw dislocations, the kinetic theory expression over¬ 
estimates the impact of GNDs on the annihilation process, whereas in microstructures with screw 
oriented GNDs, the impact of GNDs is under-estimated. 

Assuming an equi-convex microstructure, the annihilation rate of the total curvature density 
can be straightforwardly evaluated from the dislocation density annihilation rate: 


dtq' 


, = dtp' 

anil 


p' 


(B.13) 


The concomitant reduction in dislocation curvature density decreases the elongation (source) term 
vq^ in the evolution equation of the total dislocation density (461 - an effect which has important 


long-term impacts on the evolution of dislocation microstructure and may outweigh the direct effect 
of annihilation. 
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GND fraction 

Figure B.ll: Normalized annihilation rate as a function of GND fraction for screw-triggered dis¬ 
location annihilation: full line: annihilation rate for a DODF which is symmetrical around edge 
orientation (edge oriented GND), red line: polynomial approximation (1 — .2x^ — to full black 

line, dash-dotted line: result for a DODF which is symmetric around screw orientation (screw ori¬ 
ented GND), dashed line: Parabolic dependency representing the kinetic theory result for parallel 
screw dislocations. 


We have validated our model by analyzing the evolution of a system of loops with the same 
initial radius and dislocation density as in our DDD reference model. Loops are initially randomly 
distributed in a periodically continued cuboidal volume and expand with constant velocity v. Two 
loops merge if two conditions are fulfilled: (i) the contacting segments are within an angle TAy? 
from the screw orientations cp = 0,tp = tt. (ii) the spacing of the slip planes on which the loops 
are positioned is less than 1 /^. For the annihilation distance we use the value yg, = 50 nm given by 
Mughrabi et al. (19791 as typical for screw dislocation annihilation in PSB microstructures. For 


the corresponding orientation interval we use the value Acj) = ±15° given by Hussein et al. (2015). 
The results of the validation exercise shown in Fig. B.12 show an excellent agreement between 
the average evolution of dislocation density and curvature density as averaged over an ensemble of 
multiple GDD simulations and the respective predictions based upon Eqs. (B.12) and (B.13). 

To illustrate the influence of annihilation in our benchmark example we have repeated the 
simulation of the GDD^^Wodel by including the annihilation terms derived above. Parameters 
used were the same as in the validation exercise. The results shown in figure Fig. B.13 demonstrate 
that, as already estimated in the main text of the paper, for the strains and dislocation densities 
considered, annihilation does not have a very significant influence on mobile dislocation density 
evolution over a single half cycle. 
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Figure B.12: Evolution of total dislocation density and curvature density in a system of loops which 
annihilate by merging of near-screw segments. Dashed lines: averages over 200 discrete simulations, 
Full lines: CDD predictions, thin lines: evolution in individual DDD runs: for parameters see text. 



x[nm\ x[nm] 

Figure B.13: Influence of annihilation on the evolution of GND density and total dislocation density 
in CDDl; black line: result without annihilation, red line: result with annihilation. Parameters: 
A4> = 7r/12, j/a = 50 nm. 
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